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We compute the entropy density of the confined phase of QCD without quarks on the lattice to 
very high accuracy. The results are compared to the entropy density of free glueballs, where we 
include all the known glueball states below the two-particle threshold. We find that an excellent, 
parameter-free description of the entropy density between 0.7T C and T c is obtained by extending 
the spectrum with the exponential spectrum of the closed bosonic string. 
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I. INTRODUCTION 

The phase diagram of quantum chromodynamics 
(QCD) is being actively studied in heavy ion collision 
experiments as well as theoretically. A form of mat- 
ter with remarkable properties [l[ has been observed 
in the Relativistic Heavy Ion Collider (RHIC) experi- 
ments @, H, 13, Hi- It appears to be a strongly coupled 
plasma of quarks and gluons (QGP), but no consensus 
on a physical picture that accounts for both equilibrium 
and non-equilibrium properties has been reached yet. On 
the other hand, below the short interval of temperatures 
where the transition from the confined phase to the QGP 
takes place [f| 0, H, Q , it is widely believed that the most 
prominent degrees of freedom are the ordinary hadrons. 
From this point of view, the zeroth order approximation 
to the properties of the system is to treat the hadrons 
as infinitely narrow and non-interacting. We will refer to 
this approximation as the hadron resonance gas model 
(HRG). The HRG predictions were compared with lat- 
tice QCD thermodynamics data in @, , and lately they 
have been used to extrapolate certain results to zero tem- 
perature 0] ■ The HRG is also the basis of the statistical 
model currently applied to the analysis of hadron yields 
in heavy ion collisions [ll|, and recently the transport 
properties of a relativistic hadron gas have been studied 
in detail [l2j. 

Since any heavy ion reaction ends up in the low- 
temperature phase of QCD, it is important to under- 
stand its properties in detail in order to extract those of 
the high-temperature phase with minimal uncertainty. In 
this Letter we study whether the HRG model works in the 
absence of quarks, in other words in the pure SU(iV = 3) 
gauge theory, where the low-lying states are glueballs. 
There are reasons to believe that if the HRG model is 
to work at any quark content of QCD, it is in the zero- 
flavor case. Firstly, the mass gap in SU(3) gauge theory 
is very large, Mq/T c ~ 5.3. As we shall see, the thermo- 
dynamic properties up to quite close to T c are dominated 
by the states below the two-particle threshold, which 
are exactly stable. Furthermore, because of their large 
mass, neglecting their thermal width should be a good 
approximation. Secondly, the scattering amplitudes be- 



tween glueballs are parametrically 1 /TV 2 suppressed while 
those between mesons are only 1/N suppressed [13]. This 
means that the glueballs should be free to a better ap- 
proximation than the hadrons of realistic QCD. 

An additional motivation to study the thermodyna- 
mics of the confined phase of SU(3) gauge theory is that it 
is a parameter-free theory, simplifying the interpretation 
of its properties. Its spectrum is known quite accurately 
up to the two-particle threshold [3, Il5j |. By contrast, 
in full QCD calculations, lattice data calculated at pion 
masses larger than in Nature are often compared out of 
necessity to the HRG model based on the experimental 
spectrum [f| 0] ■ Finally, calculations in the pure gauge 
theory are at least two orders of magnitude faster, which 
allows us to reach a high level of control of statistical and 
systematic errors; in particular, we are able to perform 
calculations in very large volumes. 

II. LATTICE CALCULATION 

We use Monte-Carlo simulations of the Wilson action 
= ^E^,,Tr{l - P uv (x)} for SU(3) gauge the- 

ory [l6j], where is the plaquette. The lattice spacing 
is related to the bare coupling through gfi ~ 1/ log(l/aA). 
We calculate the thermal expectation value of 9 = T MM , 
the (anomalous) trace of the energy-momentum tensor 
T^y, and of #oo = ^oo ~~\®- ^ n the thermodynamic limit, 

Ts = e+p= t(0oo>T, e - 3p = (0>t - Wo- W 

Here e,p,s are respectively the energy density, pressure 
and entropy density. The operator # o = §(— E° • E a + 
B a • B a ) requires no subtraction, because its vacuum 
expectation value vanishes. The choice of of 9qo and 
9 as independent linear combinations is convenient be- 
cause they both renormalize multiplicatively. We use the 
'HYP-clover' discretization of the energy-momentum ten- 
sor introduced in [13, HH ■ The normalization of the 9oo 
operator differs from its naive value by a factor that we 
parametrize as Z(go)x(go)- The factor Z(g$) is taken 
from [h| and rests on the results of [2(3]; its accuracy 
is about one percent. The factor x(do) is obtained by 
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FIG. 1: Finite volume effects on the entropy density close to 
the deconfining temperature T c . 



Flux-Loop Mass (N c =3, N f =0) 
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FIG. 2: The mass of the temporal flux loop as calculated from 
Polyakov loop correlators, and the fit®. The Aft > 11 data 
are from the N t = 5 data from [27^7 



calibrating our discretization to the 'bare plaquette' dis- 
cretization in the deconfined phase at Nt = 6 We 
find, for 6/g% between 5.90 and 6.41, x(.9o) = 0.1306 • 
OVffo) — 0.1865 with an accuracy of half a percent. For 
the lattice beta-function that renormalizes 9, we use the 
parametrization plj of the data in [22| and the same ca- 
libration method. 

Our results for the entropy density from N t = 8 and 
Nt = 12 simulations arc shown on Fig. ([3]). The displayed 
error bars do not contain the uncertainty on the normal- 
ization factor, which is much smaller and would introduce 
correlation between the points. This factor varies by only 
7% over the displayed interval and so to a first approxi- 
mation amounts to an overall normalization of the curve. 
Our data is about five times statistically more accurate 
than that of previous thermodynamic studies [23|, [24| . 
which were primarily focused on the deconfined phase. 
Just as importantly, we kept the finite-spatial-volume ef- 
fects under good control, in particular very close to T c . 

Figure ([T]) shows the size of finite-volume effects. For 
instance, at 0.985T C the conventional choice LT = 4 leads 
to an overestimate of the entropy density by a factor 
three. The fact that the Nt = 12 data fall on the same 
smooth curve as the N t = 8 is strong evidence that dis- 
cretization errors are small. We parametrize the volume 
dependence empirically by a A+Be~ cLT curve, and use it 
to convert the N t = 12 data to LT = 8. At 0.929T C , there 
is no statistically significant difference between LT = 6 
and 8 and we do not apply any correction. It is the cor- 
rected Nt = 12 data that is then displayed on Fig. ([3]). 

In (25|, formulas for the leading finite- volume effects 
on the thermodynamic potentials were derived in terms 
of the energy gap of the theory defined on a (1/T) xLxL 
spatial hypertorus. Close to T c , this gap corresponds to 
the mass of the ground state flux loop winding around 
the cycle of length 1/T. If Ss(T, L) = s{T, oo) - s(T, L), 



the formula then reads 

-m(T)L 

Ss(T,L) = -^- r - [m 2 (T) + |T9 T m 2 (T)] . (2) 

Using the calculation of m(T) described in the next sec- 
tion, the predicted asymptotic approach to the infinite- 
volume entropy density for 0.985T C is displayed on 
Fig. ([T]). While the sign is correct, the magnitude of the 
finite- volume effects is not reproduced for LT < 8. We 
conclude that the asymptotic approach to infinite volume 
sets in for very large values of LT. Since m(T)L is only 
about 4 when LT = 6, it is not implausible that flux-loop 
states with high multiplicity dominate the finite-volume 
effects at that box size. 

Next we obtain the correlation length £(T) of the 
order parameter for the deconfining phase transition, 
the Polyakov loop. The method consists in computing 
the two-point function of zero-momentum operators, de- 
signed to have large overlaps with the ground state flux 
loop, along a spatial direction. We fit the lattice data for 
m(T) = 1/£(T) displayed on Fig. © with the formula 

( m ^ )T ) = a o~ a i(^) ~ a *{r~) ( 3 ) 

and find, either fitting a 2 or setting it to zero, 

a = 5.76(15), ai = 4.97(65), a 2 = 0.55(54) (4) 
a = 5.90(9), oi = 5.62(10), a 2 = (5) 

with in both cases a x 2 /dof of about 0.3. We remark 
that the a-i are not far from the Nambu-Goto string (28j 
values a x = = 5.02(5) [27] and a 2 = (a is the ten- 

sion of the confining string). We extract the 'Hagedorn' 
temperature, defined as in (29l | by m(Th) — 0, from the 
second fit, 

T h /T c = 1.024(3). (6) 
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This extraction amounts to assuming mean-field expo- 
nents near Th (it is not clear which universality class 
should be used [3(3]). The result is stable if the fit in- 
terval is varied, and also if 02 is fitted with ao and a\ 
constrained to the known values of (a/T 2 ) 2 and ^pyj. 

As a check on the normalization of the operators #00 
and 6, we calculate the latent heat in two different ways. 
The latent heat is the jump in energy density at T c . Since 
the pressure is continuous, we obtain it instead from the 
discontinuity in entropy density or the 'conformality mea- 
sure' e — 3p. We obtain s and e — 3p on either side of 
T c by extrapolating LT = 10 data from the confined (de- 
confined) phase towards T c . The result is 



Entropy of the confined phase (N c =3, N f =0) 



As 



1.45(5)(5), 



A(e - 3p) 



1.39(4)(5), (7) 



where the first error is statistical and the second comes 
from the uncertainty in the extrapolation (taken to be 
the difference between a linear and quadratic fit). The 
compatibility between these two estimates of Lh/T^ is 
strong evidence that we control the normalization of our 
operators. They are in good agreement with previous cal- 
culations of the latent heat performed on coarser lattices 
[27L [3l| . We have also verified more generally that the 
thermodynamic identity Td T {s/T 3 ) = (l/T 3 )<9 T (e - 3p) 
is satisfied within statistical errors. 



III. INTERPRETATION 

In infinite volume the pressure associated with a single 
non-interacting, relativistic particle species of mass M 
with n a polarization states reads 



00 -. 

P =^- 2 M>T*J2-LK 2 (nM/T) 



(8) 



where K% is a modified Bessel function. By linearity, 
the knowledge of the glueball spectrum leads to a simple 
prediction for the pressure and entropy density s = 
which is expected to become exact in the large- N limit. 
Since only the low-lying spectrum of glueballs is known, 
it is useful to consider how the density of states might be 
extended above the two-particle threshold 2Mo, where 
Mq is the mass of the lightest (scalar) glueball. The 
asymptotic closed bosonic string density of states in four 
dimensions is given by [32] 



p{M) 



( 2 ^) 3 (I* \ e M/T h 



27T h \M 



(9) 



In the string theory, the Hagedorn temperature Th is re- 
lated to the string tension, Tfi = corresponding to 
Th/T c — 1.069(5) [HI]. Below we use this value as an 
alternative to the more direct determination ([6]). 

On Fig. [31 we show the entropy contribution of the 
glueballs lying below the two-particle threshold 2Mq. 
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FIG. 3: The entropy density in units of T 3 for LT — 8. We 
applied a (modest) volume-correction to the Nt = 12 data. 



The curve is just about consistent with the smallest tem- 
perature lattice data point, but clearly fails to reproduce 
the strong increase in entropy density as T — > T c . The 
figure also illustrates that the two lowest-lying states, the 
scalar and tensor glueballs, account for about three quar- 
ters of the stable glueballs' contribution. We have used 
the continuum-extrapolated lattice spectrum [3, [34[ ■ 

Adding the Hagedorn spectrum contribution, Eq. ([9|) 
with Th given by Eq. ©, leads to the solid curve on 
Fig. [3J It describes the direct calculation of the entropy 
density surprisingly well, particularly close to T c . The 
curve tends to underestimate somewhat the entropy den- 
sity at the lower temperatures. This is likely to be a cut- 
off effect. Indeed, at fixed Nt lower temperatures corre- 
spond to a coarser lattice spacing, and the scalar glueball 
mass in physical units is known to be smaller on coarse 
lattices with the Wilson action (35|. If we use the stable 
glueball spectrum calculated at g 2 = 1 instead of the con- 
tinuum spectrum, the agreement of the non-interacting 
glueball + Hagedorn spectrum with the lattice data at 
the lower four temperatures is again excellent. This dif- 
ference provides an estimate for the size of lattice effects. 

To summarize, we have computed to high accuracy the 
entropy of the confined phase of QCD without quarks. 
The low-lying states of the theory are therefore bound 
states calledglueballs, and their spectrum is well deter- 
mined [3, [l5| . If the size N of the gauge group is in- 
creased, the interactions of the glueballs are expected to 
be suppressed [13j. To what extent the glueballs really 
are weakly interacting at N = 3 is not known precisely. 
Some evidence for the smallness of their low-energy in- 
teractions was found some time ago by looking at the 
finite- volume effects on their masses [26] . But it seems 
unlikely that glueballs well above the two-particle thresh- 
old would have a small decay width. We have neverthe- 
less compared the entropy density data to the entropy 
density of a gas of non-interacting glueballs. While re- 
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stricting the spectral sum to the stable glueballs leads to 
an underestimate by at least a factor two of the entropy 
density near T c , extending the spectral sum with an ex- 
ponential spectrum p(M) ~ exp(M/T/ l ), suggested long 
ago by Hagedorn [36|, leads to a prediction in excellent 
agreement with the lattice data for the entropy density 
(Fig. [3]). This is remarkable, since the analytic form of 
the asymptotic spectrum is completely predicted by free 
bosonic string theory, including its overall normalization 
(Eq. [5]). Therefore, since we also separately computed 
the temperature (identified with Th) where the flux loop 
mass vanishes, no parameter was fitted in the comparison 
with the thermodynamic data. By contrast, the entropy 
density is not nearly as well described if the Nambu-Goto 
value of Th is used, see Fig. 

The success of the non-interacting string density of 
states in reproducing the entropy density suggests that 
once the Hagedorn temperature has been determined di- 
rectly from the divergence of the flux-loop correlation 
length, the residual effects of interactions on the thermo- 
dynamic potentials are small. It may be that thermody- 
namic properties in general are not strongly influenced by 
interactions when a large number of states are contribut- 
ing. A well-known example is provided by the J\f = 4 
super- Yang-Mills theory, whose entropy density at very 
strong coupling is only reduced by a factor 3/4 with re- 



spect to the free theory [37|. In this interpretation, the 
main effect of interactions among glueballs on thermo- 
dynamic properties is to slightly shift the value of the 
Hagedorn temperature Th with respect to its free-string 
value. A possible mechanism is that the string tension 
that effectively determines Th is an in-medium string ten- 
sion that is ~ 8% lower than at T = 0. 

Returning to full QCD, our results lend support to the 
idea that the hadron resonance gas model can largely 
account for the thermodynamic properties of the low- 
temperature phase. Whether the open string density of 
states reproduces the entropy calculated on the lattice 
can also be tested at quark masses not necessarily as 
light as in Nature using a simple open string model [38j |. 
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